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Section  1 


INTRODUCTION 

With  the  advent  of  high  speed  computers,  complicated  and  intricate 
equations  governing  the  transfer  of  solar  radiation  in  clear  and  cloudy 
atmospheres  can  be  solved  numerically  with  high  accuracy  (see,  e.g.. 
Freeman  and  Liou,  1979).  However,  computation  time  requirements  for 
the  rigorous  solutions  to  transfer  equations  and  processes  negate  their 
usefulness  in  atmospheric  dynamic  and  climate  models  with  their  large 
number  of  grid  points  and  extended  time  integrations.  Consequently, 
new  economic  and  efficient  yet  accurate  methods  of  introducing  the 
effects  of  solar  heating  on  the  dynamics  of  the  atmosphere  must  be  sought 
to  replace  highly  theoretical  techniques. 

Simplified  methods  to  account  for  the  solar  radiation  transfer  in 
the  atmosphere  were  presented  by  Manabe  and  Moller  (1961)  and  subse¬ 
quently  by  Manabe  and  Strickler  (1964)  which  were  incorporated  in  a 
general  circulation  model  (6CM)  developed  at  the  Geophysical  Fluid 
Dynamics  Laboratory.  Moreover,  early  attempts  to  include  clouds  in 
GCM  involved  assigning  radiative  properties  to  clouds  based  upon  the 
climatological  data  such  as  those  provided  by  Haurwitz  (1948).  Katayama 
(1972)  followed  this  approach  with  a  highly  detailed  parameterization 
which  he  developed  for  the  GCM  used  at  the  University  of  California  at 
Los  Angeles.  Katayama  used  empirical  relations  to  represent  absorption 
by  water  vapor,  Rayleigh  scattering,  and  attenuation  by  dust.  His  work 
was  recently  cited  in  the  text  by  Haltiner  and  Williams  (1979).  Sasamori 


(1972)  likewise  took  this  very  simple  approach  to  the  problem  of  clouds 
when  developing  his  solar  radiation  program  which  forms  the  basis  for 
solar  radiation  input  to  the  6CM  used  at  the  National  Center  for 
Atmospheric  Research  (NCAR). 

Parameterized  solar  radiation  schemes  began  to  appear  more 
sophisticated  with  the  one  developed  by  Lacis  and  Hansen  (1974)  for 
the  GCM  developed  at  the  Goddard  Institute  for  Space  Studies  (Sommerville 
et  al.,  1974).  Their  parameterization  involved  determining  empirical 
formulae  which  accurately  fit  ozone  and  water  vapor  absorptivity  measure¬ 
ments  and  using  the  two  stream  approximation  to  determine  the  reflec¬ 
tivity  and  transmissivity  of  a  cloud  directly  from  its  optical  path. 
Multiple  scattering  effects  were  approximately  taken  into  account. 
However,  in  addition  to  the  inherent  inaccuracy  of  the  two  stream 
approximation,  their  assumption  that  all  solar  flux  was  diffuse  after 
contacting  a  cloud,  led  to  large  heating  rate  errors  below  optically 
thin  clouds  such  as  cirrus.  Additionally,  use  of  the  two  stream 
approximation  automatically  omitted  the  solar  angle  dependence  of 
cloud  transmission  and  reflection.  Recently,  Fouquart  and  Bonnel 
(1981)  presented  a  new  parameterization  for  solar  radiation  which 
is  used  in  the  GCM  at  the  Laboratory  of  Meteorology  in  Paris.  The 
results  they  show  are  generally  similar  to  those  of  Lacis  and  Hansen; 
however,  detailed  verification  of  their  method  was  not  given.  Even 
fewer  results  are  given  by  Geleyn  and  Hollingsworth  (1979)  of  their 
radiation  method  which  is  used  in  the  GCM  at  the  European  Center  for 
Medium  Range  Weather  Forecasts. 

While  each  of  the  aforementioned  schemes  appears  adequate  in  the 
context  of  computational  requirements,  verification  of  its  accuracy  is 


minimal  at  best.  Additionally,  many  simplifications  in  these  models 
appear  to  be  justified  only  by  the  results  and  have  not  been  based 
upon  physical  and  mathematical  laws  which  govern  the  radiation  field 
in  planetary  atmospheres.  Finally,  solar  radiation  interactions  with 
clouds  in  general  and  high  clouds  in  particular  have  not  been  accu¬ 
rately  treated  in  these  models. 

In  this  report,  a  new  parameterization  approach  for  the  transfer 
of  solar  radiation  in  the  atmosphere  is  developed  and  is  intended  for 
use  in  a  general  circulation  model.  Effects  of  direct  and  diffuse 
radiation  in  clear  and  cloudy  atmospheres  will  be  separately  treated, 
thus  minimizing  heating  rate  errors  which  normally  occur  below  opti¬ 
cally  thin  clouds.  Heating  rate  calculations  will  be  obtained  at 
specific  points  through  an  analytic  approach  rather  than  the  compu¬ 
tationally  more  time-consuming  and  inaccurate  finite  difference 
schemes.  Cloud  reflection  and  transmission  will  be  determined  from 
vertical  liquid  water  content  and  mean  cloud  temperature  which  affects 
the  gaseous  absorptivity  within  the  cloud.  Additionally,  a  solar 
zenith  angle  dependence  for  these  properties  will  also  be  given. 

Finally,  the  present  parameterization  model  will  allow  the  internal 
fluxes  within  cloudy  atmospheres  to  be  determined.  Ample  verification 
of  each  aspect  of  the  model  will  be  presented. 

Section  2  discusses  general  properties  of  an  absorbing  and  scat¬ 
tering  atmosphere,  justifies  the  gases  considered  in  the  model  and 
presents  scattering  characteristics  of  typical  water  and  ice  clouds. 

In  the  third  Section,  theoretical  aspects  of  computing  heating  rates 
in  a  clear  and  cloudy  atmosphere  are  shown.  The  basic  radiative 
transfer  equation  in  a  scattering  atmosphere  is  derived  and  the  discrete 
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ordinate  method  (DOM)  of  solving  it  is  discussed.  In  addition,  changes 
which  have  been  made  in  cloud  reflection  and  transmission  definitions 
given  by  Liou  (1976)  and  alterations  to  the  four  stream  DOM  program 
developed  by  Freeman  and  Liou  (1979)  are  listed.  These  char.ues  are 
required  so  that  consistent  comparisons  between  results  derived  from 
a  more  exact  radiation  program  and  those  from  the  parameterization  method 
can  be  made.  Theoretical  justifications  for  this  parameterization  model 
are  expressed  at  the  beginning  of  Section  4.  This  is  followed  by  ana¬ 
lytical  expressions  for  water  vapor,  ozone,  and  carbon  dioxide  heating 
rates  and  numerous  graphs  which  show  results  obtained  in  a  clear  atmos¬ 
phere.  Parameterization  of  solar  radiation  transfer  in  a  single  cloud 
layer  is  covered  in  Section  5.  Finally,  the  last  Section  summarizes 
the  significance  of  the  present  parameterization  method. 


The  atmosphere  is  composed  of  a  group  of  nearly  pertr''nent  gases 


and  a  group  of  gases  with  variable  concentration.  In  ac  ".ion,  it 
also  contains  various  solid  and  liquid  particles  such  as  '"osols, 

water  drops  and  ice  crystals  which  are  highly  variable  i  ;e  and 
time.  Table  1  lists  the  concentration  of  the  permanent  and  variable 
gases  in  the  earth's  atmosphere  as  shown  by  Liou  (1980). 

The  mixing  ratio  of  each  permanent  gas  is  constant  from  the 
earth's  surface  up  to  an  altitude  of  about  60  km;  however,  it  should 
be  recognized  that  the  concentration  of  carbon  dioxide  actually  varies 
as  a  result  of  the  combustion  of  fossil  fuels  and  photosynthesis  proc¬ 
esses.  The  concentration  of  water  vapor  and  ozone  are  highly  variable 
in  space  and  time.  Water  vapor  density  generally  decreases  with  height 
in  the  troposphere  while  ozone  has  a  maximum  concentration  in  the  15- 
25  km  region. 

2.2  Gaseous  Absorption 

Absorption  spectra  due  to  electronic  transitions  of  molecular  and 
atomic  oxygen  and  nitrogen,  and  ozone  occur  in  the  ultraviolet  (UV) 
region  of  the  solar  spectrum  while  those  due  to  vibrational  and 


rotational  transitions  of  triatomic  molecules  such  as  H^O,  CO2.  and 
ozone  lie  in  the  infrared  (IR)  region.  Figure  1  depicts  the  solar 
spectrum  measured  at  the  top  of  the  atmosphere  (Thekaekara,  1974), 
the  earth's  surface,  and  the  areas  of  atmospheric  absorption.  Absorp¬ 
tion  of  UV  radiation  occurs  primarily  in  the  stratosphere  and  mesosphere 
while  energy  in  the  longer  wavelengths  is  absorbed  in  the  troposphere. 
Of  the  gases  in  the  upper  atmosphere,  only  ozone  strongly  absorbs 

O 

at  wavelengths  above  2000  A  where  significant  solar  energy  is  found. 

O 

The  ozone  absorption  spectrum  begins  at  2000  A  and  continues  to  11800 

o  o 

A.  The  area  of  strongest  absorption  occurs  between  2000-3000  A  in 
the  Hartley  band  complex.  Another  area  of  strong  ozone  absorption 

O 

is  found  in  the  Huggins  band  complex  located  in  the  3000-3600  A  region. 
These  bands  are  characterized  by  highly  variable  absorption  coeffi- 

O 

cients.  The  much  weaker  Chappius  bands  absorb  in  the  4400-11800  A 
band  interval. 

Water  vapor  is  the  most  important  absorber  in  the  near  infrared 
region  of  the  solar  spectrum  and  it  almost  exclusively  is  responsible 
for  atmospheric  heating  below  10  kilometers.  This  absorption  occurs 
in  six  bands.  They  are  the  overtone  and  combination  bands  centered 
at  0.94,  1.1,  1.38,  and  1.87  pm,  the  2.7,  and  3.2  pm  bands.  Although 
the  2.7  pm  band  is  the  most  important,  significant  absorption  also 
occurs  in  the  3.2  pm  band  and  the  overtone  and  combination  bands. 

One  additional  water  vapor  absorption  band  is  located  at  6.3  pm  where 

a  negligible  solar  flux  is  observed. 

Carbon  dioxide  exhibits  a  number  of  weak  absorption  bands  in  the 
near  IR  region  of  the  solar  spectrum.  Those  centered  at  5.1,  4.8, 

4.3,  2.0,  1.6,  and  1.4  pm  are  very  weak  and  can  be  ignored.  The  2.7 
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vn  band  is  somewhat  stronger.  Carbon  dioxide  absorption  is  compli¬ 
cated  in  this  band  since  it  overlaps  with  the  2.7  water  vapor  band. 


Oxygen  absorbs  very  strongly,  as  does  nitrogen,  at  wavelengths 

O 

smaller  than  2000  A  but  Park  and  London  (1974)  have  demonstrated  that 
it  is  of  minor  importance  in  heating  rate  calculations  below  75  kilo¬ 
meters.  Since  the  area  of  interest  in  this  study  extends  to  only  60 
km,  this  research  will  concentrate  on  gaseous  absorption  by  water 
vapor,  ozone  and  carbon  dioxide. 


2.2.1  Absorption  by  Ozone 

In  the  15-60  km  region,  atmospheric  heating  results  primarily 
from  the  absorption  of  solar  radiation  by  ozone.  Ozone  absorption 
occurs  in  two  complex  band  systems  as  shown  in  Figs.  2  and  3.  In 

o 

the  Hartley  and  Huggins  bands  (2000-3600  A),  ozone  absorbs  very 
strongly  while  it  absorbs  much  less  strongly  in  the  Chappius  bands 
(4400-11800  A). 

Ozone  absorption  was  modeled  by  an  absorbing  layer  overlaying  a 
reflecting  layer  as  described  by  Lacis  and  Hansen  (1974).  The  spectral 
absorptivity  for  the  solar  flux  having  a  cosine  of  solar  zenith  angle 
Pq  and  corresponding  to  an  ozone  path  length  u^  (cm  NTP)  is 


VI 


=  1  -  e 


-kjUj/uo 


(1) 


where  k^  is  the  monochromatic  ozone  absorption  coefficient  after  Inn 
and  Tanaka  (1953)  and  Vigroux  (1953)  as  compiled  by  Howard  (1961). 
These  data  are  listed  in  Table  2.  Absorption  coefficients  in  the 
Hartley  and  Huggins  band  complex  have  a  temperature  dependence 
(data  listed  correspond  to  -44°  C)  while  those  in  the  Chappius  band 


Fig.  3.  Absorption  spectrum  for  ozone  in  the  Chappius  band  complex  from  4000  to  11800  A. 


2.  Monochromatic  Absorption  Coefficient: 
Ozone. 


X(vim) 

.20 

.21 


.22 


.23 

.24 

.25 

.26 

.27 


.28 


.30 


.32 


.34 


.36 


.45 


.50 


.55 

.60 


.65 


.70 


.80 


A  (cm”^) 
v'  ' 

8.61 
14.7 
48.4 
122.0 
216.0 
299.0 
292.0 
210.0 
106.0 
10.1 
.898 
6.40x10"^ 
1.80x10"^ 
3.50x10"^ 
3.45x10'^ 
9.20x10"^ 
.132 

6.20x10“^ 

2.30x10"^ 

1.00x10"^ 
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complex  have  practically  no  temperature  dependence.  Absorption  coeffi- 

O 

cients  for  wavelengths  greater  than  8000  A  are  very  small  and  are 
neglected  in  ozone  heating  calculations. 

A  few  words  are  in  order  about  the  units  used  to  express  ozone 
density,  since  it  can  be  quite  confusing.  The  familiar  cgs  units  of 

_3 

gm  cm  are  seldom  used  in  ozone  calculations  but  an  older  unit  which 

-3  -1 

is  frequently  used  is  the  10  cm  NTP  km  .  The  conversion  is 

03(9  cm"^)  =  21.4  P3  {10‘^  cm  NTP  km"^)  .  (2) 


2.2.2  Water  Vapor  Absorption 

Howard  et  al.  (1956)  published  results  of  extensive  laboratory 
experiments  which  measured  the  total  absorption  of  the  0.94,  1.1,  1.38, 
1.87,  2.7,  and  3.2  ym  water  vapor  bands  for  various  path  lengths  and 
under  different  atmospheric  pressures.  They  derived  two  formulas  for 
this  total  absorption.  For  small  values  of  total  absorption 


A  =  I  A^dv  =  cuj'^^  (Pg  +  e)*^  ,  A  >  A^ 


(3) 


and  for  large  values  of  total  absorption 


A  =  A,dv  =  C  +  0  log  u,  +  K  log  (P.  +  e)  ,  A  >  A^  (4) 

V  X  Q  c 

4 

where  v  is  the  wavenumber  (cm”^),  A  the  band  area  (cm"^),  A^  the  frac¬ 
tional  absorption  within  the  band  at  v,  u^  the  water  vapor  path  length 
.0 

(g  cm  ),  e  the  partial  pressure  of  the  absorbing  gas  (mm  Hg),  P^  the 
partial  pressure  of  the  non-absorbing  gases  (mm  Hg),  A^  the  critical 
band  area  above  which  the  strong  band  equation  [Eq.  (2)]  applies, 
and  c,  k,  C,  D,  and  K  are  empirical  constants  listed  in  Table  3. 
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Since  these  two  formulas  were  not  continuous  when  A  =  A  ,  Liou  and 

w 

Sasamori  (1975)  derived  a  single  formula  to  approximate  the  mean  band 
absorptivity,  A-,  for  both  weak  and  strong  absorptivity  in  the  follow¬ 
ing  form: 

A-  =  A/Av  =  -^  [C  +  D  logjQ  (x  +  Xq)]  ,  (5) 

where 

X  =  Uj  .  (6) 

Xq  =  .  (7) 

and 

P  *  f  ^  P(u)  du  /  du  .  (8) 

Jo  Jo 

The  reduced  pressure  accounts  for  the  pressure  dependence  of  absorp¬ 
tion  in  a  non-homogeneous  atmosphere.  It  is  an  empirical  method  and 
places  all  the  absorbing  matter  along  a  pressure  gradient  at  one 
pressure, 

2.2.3  Absorption  by  Carbon  Dioxide 


As  previously  indicated  in  Table  1,  CO^  was  assumed  to  be  a 
permanent  constituent  in  the  atmosphere  with  a  mixing  ratio  of  330 
parts-per-million  (ppm)  by  volume.  It  absorbs  solar  radiation  most 
strongly  in  the  2.7  vm  band  region  where  it  overlaps  with  water  vapor. 
Other  minor  CO2  absorption  bands  occur  at  1.4,  1.6,  2.0,  4.3,  4.8,  and 
5.2  v<n,  but  were  neglected  because  absorption  is  very  weak  in  compar¬ 
ison  with  water  vapor  absorption. 


Howard  et  a1 .  (1956)  also  measured  the  total  absorptivity  of  each 
CO2  band  and  developed  formulas  to  approximate  it.  Similarly,  Liou 
and  Sasamori  (1975)  combined  these  formulae  into  Eq.  (5).  The  per¬ 
tinent  constants  are  listed  in  Table  4  for  CU^.  Absj'  O?'^  path  h  ,vgih 
units  for  COg  are  in  cm-atm. 

In  another  of  the  series  of  papers  published  on  gaseous  absorp¬ 
tion  by  Howard  et  al .  (1956),  they  theorized  the  accepted  method  for 
computing  the  total  absorptivity  due  to  water  vapor  and  carbon  dioxide 
overlapping  in  the  same  spectral  region.  The  total  gaseous  absorption 
in  the  2.7  pm  water  vapor  and  carbon  dioxide  overlap  band  was  given  as 


A^(H20  +  CO2)  dv  =  A^(H20)  dv  +  [  6  A  (CO2)  dv 


V  =  2.7  pm 


where  6  is  a  weighting  factor  which  ranges  between  0  and  1,  depending 
on  the  total  water  vapor  absorptivity  in  the  2.7  pm  band.  Its  values 
are  given  in  Fig.  4.  The  values  of  6  were  calculated  using  one  of  the 
following  equations  after  assuming  Fig.  4  could  be  described  accu¬ 
rately  by  two  straight  lines.  If  the  total  water  vapor  absorption, 
Ag(ui/Po)  Avg,  in  the  band  is  greater  than  or  equal  to  300  cm"^, 

^  -  §55  '"'5  " 


Otherwise 
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Fig.  4  Carbon  dioxide  weight  factor  as  a  function  of  total 
water  vapor  absorption  in  2.7  pm  band. 
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Section  3 


RADIATIVE  TRANSFER  IN  THE  ATMOSPHERE 


3.1  General 

The  exact  solution  to  the  radiative  transfer  equation  can  be 
computed  with  high  accuracy  if  computational  time  constraints  are 
not  important.  This  section  discusses  the  derivation  of  the  radia¬ 
tive  transfer  equation  (RTE)  for  diffuse  intensity  in  a  plane-parallel 
atmosphere,  the  discrete-ordinate  method  (DOM)  of  approximating  the 
exact  solution  to  the  RTE,  the  four  stream  DOM  analytical  solution 
and  modifications  made  to  the  University  of  Utah's  four  stream  DOM 
program  used  for  comparison  with  the  parameterized  model. 

The  concepts  discussed  in  the  derivation  of  the  RTE  for  diffuse 
intensity  in  a  plane-parallel  atmosphere  follow  closely  the  description 
given  in  Liou  (1980). 

3.2  Derivation  gf  the  Radiative  Transfer  Equation  for  Diffuse 

Intensity  in  a  Plane-Parallel  Atmosphere 

In  a  plane-parallel  atmosphere,  the  differential  diffuse  inten¬ 
sity  in  a  directional  stream  of  radiation  defined  by  a  solid  angle 
is  reduced  by  single  scattering  of  energy  out  of  the  direction  ^ 
and  absorption  by  gases  and  particles.  This  is  expressed  by 

dl(z,i^)  =  -  Og  ^  I(z,j^)  dz  ,  (12) 

where  dz  is  the  differential  thickness  of  a  layer,  the  mean  extinction 
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cross  section  of  a  particle,  N  the  total  number  of  particles  per  volume, 
and  ^  the  directional  element  of  solid  angle  that  represents  the  pencil 
of  radiation  (see  Fig.  5).  is  described  by  azimuth,  (j),  and  zenith, 
n,  angles.  Meanwhile,  the  differential  diffuse  intensity  in  tie  i: 
ection  ^  may  be  increased  by  multiple  scattering  of  radiation  in  che 
direction  n' .  This  is  given  by 


dl(z,^) 


0$  N 
cose 


HzM‘) 


P(^,^' ) 

4it 


dfi'  , 


(13) 


where  is  the  mean  scattering  cross  section  and  the  phase  function, 
P(S^,S1'),  introduces  the  appropriate  radiation  stream  from  to 
Thus  integration  over  all  solid  angles  gives  all  possible  contribu¬ 
tions  of  multiply  scattered  energy  from  SI'  to  Finally,  the 
differential  diffuse  intensity  in  direction  n  may  be  increased  by 
single  scattering  of  the  direct  solar  radiation  whose  direction  is 
represented  by -^Q,  where  the  minus  sign  denotes  the  downward  direc¬ 
tion  of  the  direct  solar  beam.  If  the  direct  solar  flux  arriving  at 
level  z  is  F(z),  then  this  effect  is  given  by 

o  Ndz  P(^,-fin) 


where 


F(z)  =  ttFq  exp  ^ 


cosG, 


0^(z')  N(z')  dz'> 


(15) 


and  ttFq  is  the  flux  of  solar  energy  arriving  at  the  top  of  the  atmos¬ 
phere  at  angle  Og  to  the  local  zenith.  Upon  combining  these  three 
terms  and  introducing  the  optical  depth,  t,  we  get 
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o,b  Diffuse 
c  Dired 


Transfer  of  solar  radiation  in  plane-parallel  layers, 
a:  attenuation  by  extinction,  b;  multiple  scattering, 
c:  single  scattering  of  the  unscattered  solar  flux. 


(16 


where 


4^  ttFq  ® 


* 

=  a  (z*)  N(z')  dz'  , 
J7 


dx  =  -  Og  N(z)  dz  • 


The  single  scattering  albedo  5  is  defined  as  the  ratio  of  scattering 
cross  section  to  extinction  cross-section  o^.  In  this  report,  we 
describe  for  documentation  purposes  the  discrete-ordinates  method 
for  the  solution  of  the  radiative  transfer  equation  governing  the 
inhomogeneous  gaseous  atmosphere. 


3.3.  The  Discrete-Ordinates  Method 

3.3.1  Intensity  in  a  Non-Homogeneous  Atmosphere 


Since  only  fluxes  are  needed  in  the  calculation  of  solar  heating 
rates,  Eq.  (16)  may  be  written  for  the  azimuthally  independent  dif¬ 
fuse  intensity  of  solar  radiation  in  a  plane-parallel  homogeneous 
layer,  as 
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(19) 


£  +1 

y  ^  =  i^(t.p)  -  I  I  j  I^(t.p')  dp' 


Liou  (1973)  has  shown  the  analytical  solutions  for  Eq.  (19)  based 
upon  the  discrete  ordinate  method  for  radiative  transfer  are  given  by 


I  (t.p^) 


I  ^  z^(Pt)  exp  (--[/Po) 

^  t  =  1 . N. 

i  =  -n,  .  .  ,  ,  n  .  (20) 


In  Eq.  (20)  ^  denotes  j  from  -n  to  n  (n  +  0),  and  the  eigenfunctions 

j 

are  derived  from  the  associated  homogeneous  system  as 


4>^(p-)  “  — ^'T  y  C  (1^^)  P^(p*)  » 

J  '  1+p.k^  m=0  "•  ">  J  "> 

1  J 


where  M  denotes  the  number  of  terms  in  the  Legendre  polynomial,  p^, 
expansion,  is  a  set  of  M  +  1  constants,  and  is  a  constant  of 


proportionality  and  can  be  determined  by 


3m+l  -  u*'  „ 

- 01  ^  -  -'"t  €  1  .  m  =  0.1, 

k^mn) 


.  .,  M-1  . 


By  taking  Cq  =  1»  ^ll  values  of  ^  can  be  evaluated.  The  eigenvalues, 
k  ,  may  also  be  determined  from 


f(k*')  =  1  -  I  I  =  0  . 
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The  last  term  in  Eq.  (20)  represents  a  particular  solution  where 


fl  1  ^  VPa/  ”  '  •!  9  f  1  ^ 

*^0  -  - - Jo 


and  the  H  function  is  given  by 


8  1  "  "  £ 

"  u  L  —  "iT  ^  (w+Pi)  Tt  (1  +  k.w) 
UlP2*--Pn  i=i  ^  j=l  J 


The  solution  expressed  in  Eq.  (20)  is  valid  only  for  nonconservatf ve 
scattering.  For  conservative  scattering  a  slightly  different  form  is 
needed  but  neglected  here  for  simplicity. 

The  Lj  in  Eq.  (20)  are  a  set  of  N  x  2n  coefficients  to  be  deter- 
minded  from  boundary  and  continuity  equations  in  a  non -homogeneous 
atmosphere.  The  atmosphere  is  inhomogeneous  in  that  the  single¬ 
scattering  albedo  and  phase  function  are  generally  nonlinear  functions 

of  the  optical  depth.  To  account  for  this  inhomogeneity,  it  is 
customary  to  divide  the  atmosphere  into  a  number  of  layers  such  that 

each  layer  may  be  considered  to  be  optically  homogeneous. 

The  boundary  conditions  which  describe  the  diffuse  intensity  for 
solar  radiation  at  the  top  and  bottom  of  the  atmosphere  are,  respec¬ 
tively. 


l\o,p.)  =0  i  =  1,2,  .  .  .,  n  ,  (26) 

I  (t  ,+u^)  -  [F  (t  )  +  Pq  IT  Fq  exp  (-T  /Pq)]  , 

i  =  1,  2,  .  .  .,  n  ,.  (27) 
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where 


F  (t  )  -  271  I  I  (7  ,-p.)  a.  p. 

1=1  111 


with  being  the  Gaussian  quadrature  weights  and  assuming  an  Iso¬ 
tropic  reflecting  surface.  These  conditions  state  that  there  is  no  down¬ 
ward  (-p^)  diffuse  intensity  at  the  top  of  the  atmosphere,  while  at 
the  bottom  of  the  atmosphere,  the  upward  (+p^)  diffuse  intensity  is 
the  product  of  the  surface  albedo  (r^)  and  the  total  flux  reaching 
the  surface  t  . 

The  continuity  condition  for  atmospheric  layers  is 


,1,  I  X  .1+1,  t.  X 

I  (t  ,  ,p.)  =  I  (i  ,p^)  ,  1  =  -n, 


. ,  n  , 


£  =  1,  .  .  .  N-1  .  (29] 

This  condition  requires  no  energy  loss,  in  either  the  upward  or  down¬ 
ward  radiation  streams,  at  the  interface  separating  each  homogeneous 
layer.  Note  that  t  represents  the  optical  depth  from  the  top  of  the 
atmosphere  to  the  bottom  of  the  I  layer. 

Substituting  Eq.  (20)  into  Eqs.  (26)  -  (28)  produces  a  set  of 

£ 

equations  for  determining  the  L.  constants.  They  are 

I  =  -  2*'(-P^)  ,  i  =  1.2,  .  .  .,  n  ,  (30) 


I[L;y;(Pi).Lj^l  (p,)] 


n  (u^)  , 


—n*  •  >  • 

=  1,2,  ...  ,  N-1 


V  |N  oN,  »  N,^  X 

I  Lj  6j(m,)  =  -e  (‘u,)  .  1 


.  ,  n 


where 


irj(|i,)  -  ♦](!*, )  exp  (-kjx‘) 


'  [♦3(»u,)  -  2  Tg  J  ♦Jc-p,)  p,] 

J  ~  ^ 

X  exp  (-kjT^)  . 


_  r_4 


[z  (y^)  -  Z  (p.)]  exp  (-tVpq) 


e'^{+M^)  =  [2^(+y^)  -  2  Tg  z"(-yj) 


■  ^  ^0  *"0^  ®*P  (-tVo^- 


Hence,  there  are  N  x  2  n  equations  for  the  determination  of  N  x  2n 

£ 

unknown  constants  L..  In  terms  of  matrix  operations,  this  system  of 

V 

equations  becomes 


♦  L  =  X  . 


where 
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where  blank  spaces  denote  zero  elements.  With  this  Information,  the 

a 

coefficients,  L.,  may  be  determined  simultaneously  by  employing  any 

jj 

matrix  Inversion  method.  Substituting  the  L.  back  Into  Eq.  (20)  pro- 

J 

duce.'  the  intensity  field  within  each  layer. 

The  modified  radiative  transfer  program  now  allows  21  levels  and 
20  layers  in  the  atmospheric  model.  It  can  use  any  of  the  five  clima¬ 
tological  atmospheric  profiles  compiled  by  McClatchey  et  al .  (1972). 
After  the  optical  depths  are  calculated  for  each  sub-band  at  each  of 
the  21  levels,  solutions  defined  by  Eq.  (20)  are  computed  for  each  level 
In  the  atmosphere  and  for  each  of  the  four  discrete  streams  of  radia¬ 
tion.  There  is  a  separate  set  of  solutions  for  each  sub-band;  conse¬ 
quently,  the  L's  shown  In  Eq.  (20)  must  be  determined  42  times,  24  times 
for  water  vapor  and  18  times  for  ozone.  Each  of  these  42  sets  of  L's 
are  found  by  inverting  a  80x80  matrix.  After  computing  the  L's,  the 
upward  and  downward  Intensities  of  the  radiation  stream  at  a  given  level 
can  then  be  calculated. 

3.3.2  Fluxes  and  Heating  Rates 

Once  the  intensities  have  been  determined  for  each  layer,  the  sub¬ 
band  upward  and  downward  fluxes  are  computed  for  each  layer  (see  Figure 
6).  The  fluxes  are  given  by 

F^t)  =  2-n  \  a.  y.  I(t,m.)  ,  (42) 

1=1  ^  ^  ^ 

for  upward  fluxes,  and 

F  (t)  =  -  2t[  \  a.  y.  I(T,-y.)  -  S(T,yQ)  , 


(43) 


.  Configuration  of  solar  radiation  flux  transfer  in  cloudy  atmospheres 


for*  downward  fluxes,  where  S(i,Uq)  is  the  direct  solar  flux  arriving 
at  optical  depth  t/Mq  and  is  given  by 


S(i,ijq)  =  Uq  Fq  exp  (-i/Uq)  . 


The  net  flux  for  a  given  layer  may  be  computed  from 


F(c)  =  F^t)  +  F^t)  , 


and  the  sub-band  heating  rate  for  a  layer  of  air  resulting  from  the 
absorption  by  atmospheric  gases  and  cloud  particles  for  any  given 
band  can  be  calculated  from 


E]  - 

3Z 


1 


where  p.  is  the  mean  density  of  air  in  layer  I  and  is  the  specific 
a  p 

heat  of  air  for  constant  pressure.  The  sub-band  heating  rates  are 
summed  together  and  weighted  according  to  the  solar  flux  in  each  sub¬ 
band  to  produce  the  total  heating  rate  for  a  layer. 


3.3.3  Cloud  Optical  Properties 


Clouds  regularly  occupy  about  50%  of  the  sky  on  a  global  scale 
and  are  the  single  most  important  regulators  of  the  earth's  radiation 
balance.  They  absorb  and  scatter  incoming  solar  radiation  at  all 
wavelengths.  However,  introduction  of  clouds  into  an  atmospheric 
profile  vastly  complicates  the  radiative  calculations  because  the 
scattering  source  function  in  the  transfer  equation  involves  Mie 
scatter  which  is  a  complex  function  of  cloud  microphysics. 
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To  obtain  the  reflection,  absorption,  and  transmission  of  the 
cloud  for  the  entire  solar  spectrum,  proper  summation  over  the  fluxes 
in  each  spectral  region  is  required. 

Consider  a  cloudy  atmosphere  containing  one  cloud  layer  with  top 
and  base  heights  denoted  as  and  as  shown  in  Figure  6  and  assume 

the  upward  and  downward  fluxes  are  known  at  each  level.  The  cloud 
reflection  is  defined  as  the  ratio  of  the  reflected  flux  density  to  the 
incident  solar  flux  density  normal  to  the  cloud  top.  Thus 


A.L. 

where  f^  denotes  the  fractional  solar  flux  in  the  i  spectral  band. 

The  cloud  transmission,  t^,  is  defined  as  the  ratio  of  the  downward 
flux  density  at  the  cloud  base  to  the  incident  solar  flux  normal  to  the 
cloud  top  and  is  given  by 

‘c  '  I  i  ^  •  <'>81 


Finally,  the  cloud  absorption,  is  defined  as  the  ratio  of  the  net 
flux  density  divergence  of  the  cloud  layer  to  the  incident  solar  flux 
normal  to  the  cloud  top.  Thus, 


Heating  rates  for  cloud  layers  are  computed  as  previously  des¬ 
cribed  for  clear  layers,  once  the  upward  and  downward  fluxes  are 
known  at  the  cloud  top  and  cloud  base. 
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It  should  be  noted  that  these  cloud  radiative  properties  defini¬ 
tions  are  not  the  same  as  those  previously  defined  by  Liou  (1976)  and 
Liou  and  Wittman  (1979)  which  defined  broadband  cloud  radiative  proper- 
ti  v.  Miing  inaccurate  ‘i^^ummations  over  wavelength. 

3.3.4.  The  Four  Stream  Solution 

For  discrete  streams  of  16,  (n  =  8),  the  discrete  ordinate  method 
(DOM)  produces  results  which  are  generally  accurate  up  to  3-4  digits 
when  compared  to  doubling  computations  which  have  an  accuracy  of  about 
five  decimal  points  according  to  van  de  Hulst  and  Grossman  (1968).  The 
discrete  ordinate  method  has  been  found  to  be  quite  satisfactory  because 
it  allows  the  solutions  to  Eq.  (20)  to  be  explicitly  derived  and  can 
easily  take  into  account  various  distribution,  thickness,  and  cloud 
types  making  the  scheme  well  suited  to  studying  climatic  effects  of  the 
clouds. 

Liou  (1973)  found  the  accuracy  of  results  produced  by  the  four 
stream  (n  =  2)  discrete  ordinate  method  to  be  within  about  1%  of  those 
generated  by  16  streams  for  flux  calculations;  thus,  the  four  stream 
DOM  model  developed  by  Freeman  and  Liou  (1979)  was  considered  a  more 
exact  solution  for  comparing  with  the  parameterized  model.  However, 
since  the  previous  solar  transfer  program  developed  at  the  University  of 
Utah  had  a  different  emphasis  than  the  present  study,  several  modifica¬ 
tions  were  necessary  before  it  could  be  properly  used  for  comparison 
purposes.  These  modifications  include: 

(1)  In  Freeman  and  Liou's  (1979)  work  the  objective  of  the 
research  was  not  centered  upon  ozone  heating  rates.  Consequently,  the 
radiation  model  extended  up  to  only  30  km  and  contained  just  two  absorption 
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coefficients  for  ozone;  one  in  the  Hartley-Huggins  band  complex,  another 
for  the  Chappius  band  complex.  This  model  was  modified  to  allow  12 
absorption  coefficients  in  the  former  and  six  in  the  latter  ozone  band 
complexes  while  the  height  was  extended  to  60  km.  The  ozone  absorption 
wavelength  interval,  mean  absorption  coefficient,  and  fractional  solar 
flux  in  each  of  these  18  bands  are  listed  in  Table  5.  The  mean  ozone 
absorption  coefficient  in  each  wavelength  interval  in  Table  5  was 
calculated  by  weighting  the  mean  absorption  coefficient  in  each  0.05  pm 
sub-interval  with  the  fractional  solar  flux  in  the  sub-interval,  then 
dividing  by  the  total  solar  flux  in  the  interval.  Note  that  Table 
4  only  depicts  the  monochromatic  absorption  coefficient  at  the  wave¬ 
lengths  shown. 

(2)  The  previous  solar  radiation  program  developed  at  the  University 
of  Utah  only  computes  heating  rates  in  a  clear  atmosphere  or  an  atmos¬ 
phere  containing  one  cloud.  Considerable  programming  was  required  to 
allow  multiple  cloud  layers  to  be  introduced  into  the  model  atmosphere. 

The  present  program  will  now  allow  up  to  four  cloud  layers  in  an  atmos¬ 
phere. 

(3,  The  six  water  vapor  bands  were  each  subdivided  into  five 
sub-bands  in  which  the  absorption  coefficient  was  assumed  constant. 

Thus,  the  optical  depth  was  calculated  for  each  wavelength  interval 

at  each  level  in  the  atmosphere  using  Eq.  (17).  However,  water 

vapor  absorption  coefficients  in  an  earlier  version  of  the  solar 

radiation  program  were  recomputed.  This  because  wavelength 

intervals  larger  than  those  derived  by  Liou  ui.d  Sasamori  (1975)  were 

used  in  the  present  program  which  was  intended  to  allow  absorption 

and  scattering  processes  to  be  considered  in  the  transfer  calculation.  For 
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Tdble  5.  Ozone  Wavelength  Interval  with  Associated  Mean  Absorption 


Coefficient  and  Fractional  Solar  Flux 


Wavelength 

Interval  (pm) 

.?0-.2\ 

Absorption  ^ 

Coefficient  (cm'^] 

9.8 

Fractional 

Solar  Flux 

1.24x10"^ 

.21-. 22 

27 

2.97x10"^ 

.22-. 23 

75 

4.599x10“^ 

.23-. 24 

164 

4.59x10'^ 

.24-. 25 

254 

5.14x10“^ 

.25-. 26 

290 

7.55x10"^ 

.26-. 27 

241 

1.35x10"^ 

.27-. 28 

145 

1.59x10’^ 

.28-. 30 

33.7 

6.46x10'^ 

.30-. 32 

2.8 

1.01x10'^ 

.32-. 34 

.16 

1.50x10'^ 

.34-. 35 

.014 

7.95x10"^ 

.45-. 50 

.011 

7.46x10'^ 

.50-. 55 

.055 

6.78x10"^ 

. 55- . 60 

.11 

6.30x10'^ 

.60-. 65 

.09 

5.87x10"^ 

.65-. 70 

.038 

5.33x10"^ 

.70-. 80 

.015 

9.14x10“^ 

_ _ _ 

34 


comparisons  between  results  computed  from  parameterization  and  more 
exact  transfer  programs  we  have  derived  a  new  set  of  equivalent  absorp¬ 
tion  coefficients,  k^,  and  weighting  factors,  w^,  for  each  sub-band 
in  the  six  water  vapor  bands  and  the  larger  band  frequency  interval 
which  appear  in  Table  6. 

(4)  Finally,  the“def inition  of  cloud  absorptivity,  reflectivity, 
and  transmissivity  used  by  Liou  and  Wittman  (1979)  was  modified  to 
coincide  more  accurately  with  the  broadband  parameterization  approach. 
The  broadband  definitions  for  cloud  radiative  properties  are  given  in 
Eqs.  (47)  -  (49). 
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Section  4 


THE  BROADBAND  PARAMETERIZATION  APPROACH 
IN  A  CLEAR  ATMOSPHERE 

4.1  General 

In  the  previous  section,  the  four  stream  DOM  solution  for  the  RTE 
was  described.  This  type  of  solution  is  quite  time  consuming  and 
expensive  since  heating  rates  must  be  computed  for  each  wavelength 
interval  then  integrated  across  the  frequency  spectrum  to  produce  the 
total  beating  rate  for  a  layer.  The  broadband  approach  simplifies 
this  procedure  by  integrating  the  absorption  coefficient  over  wave¬ 
length,  producing  an  absorptivity  which  is  wavelength  independent,  so 
the  RTE  need  only  be  solved  once. 


4.2  Gaseous  Absorption  in  Clear  Skies 

In  clear  skies,  atmospheric  heating  results  basically  from  gas¬ 
eous  absorption,  thus  the  radiative  transfer  equation  simplifies  to 


dl^(u) 

^0  k^du  "  "  ’ 


(50) 


where  du  is  the  path  length  of  the  absorbing  gas.  The  solution  to  this 
equation,  known  as  Beer's  Law,  is 


/  \ 

,  fU 

I^  (u)  =  1^(0)  exp  -  J  k^  du  /  , 
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which,  after  integrating  over  solid  angles  gives 


F^^(u)  =  F^(0)  exp 


-  fn  "x 


(52) 


where  I;^(0)  and  F^(0)  denote  the  monochromatic  intensity  and  flux  density 
respectively,  at  u  =  0. 

The  downward  flux  density  along  the  solar  beam  is  found  by  inte¬ 
grating  the  monochromatic  flux  density  across  wavelength  such  that 


f/(0)  exp 


►  ’ 

-  j/x 


dA  . 
(53) 


The  monochromatic  absorptivity 


in  a  non-scattering  medium  is  defined  by 


A^(u)  =  1  -  exp 


(54) 


Thus,  Eq.  (53)  becomes 


F^u)  =  1^  f/(0)  [1  -  A^(u/yQ)]  dA  .  (55) 

Defining  the  solar  constant  such  that 

.oo 

f/(0)  dA  ,  (56) 

0 

then  the  downward  flux  density  perpendicular  to  the  layer  of  strati¬ 
fication  is  given  by 

F^u)  =  p.  S.  [1  -  A(u/y.)3 


(57) 


The  heating  rate  at  level  z(u)  is  given  by 


il  -  -JL-  iE  _p_  H 

■  ■  CpPa  ^  ’ 


where  p.  is  the  air  density  at  level  z,  p  and  u  apply  to  the  absorb- 

O 

ing  gas. 

Thus,  the  broadband  concept  requires  a  gaseous  absorptivity  which 
is  wavelength  independent  and  applies  to  the  entire  solar  spectrum. 

This  may  be  achieved  by  noting  that 

A(Uj/yQ)  =  I 

where  f^.  is  the  fractional  solar  flux  for  the  i*"”  wavelength  interval 
and  is  the  monochromatic  gaseous  absorptivity  given  in  Section  2 
for  each  gas  type.  Subscript  j  denotes  the  type  of  absorbing  gas  where 
1,  2,  and  3  refer  to  water  vapor,  carbon  dioxide,  and  ozone,  respectively. 


4.2.1  Ozone  Absorption  and  Heating  Rate  Profiles 

The  broadband  ozone  absorptivity  is  given  by  Eq.  (59)  where 
i  =  1,  .  .  .  ,18  and  the  total  ozone  absorptivity  is  given  by 


A(u3/uq)  =  [1  -  exp  (-^^03/^0)^ 


where  is  the  mean  absorptivity  and  f.  is  the  fractional  flux  in  the 
i^^  spectral  interval  as  given  in  Table  6. 

In  Fig.  7,  the  ozone  broadband  absorptivity  computed  by  Eq.  (59) 
for  path  lengths  of  0.0001  to  10  cm-NTP  are  compared  with  other  pub¬ 
lished  data.  The  broadband  absorptivity  compares  very  favorably  with 
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the  line-by-line  calculations  taken  from  Fouquart  and  Bonnel  (1980)  and 
also  is  in  good  agreement  with  that  of  Lacis  and  Hansen  (1974),  although 
slightly  larger  at  very  small  path  lengths. 


The  ozone  heating  rate  may  be  determined  in  using  Eq.  (58).  The 
rate  of  change  of  downward  solar  flux  is  readily  calculated  by  substi¬ 
tuting  Eq.  (55)  into  Eq.  (58).  This  gives 


18  _ 

fj  k,  exp  (-kjUj/Pp)  . 


(61) 


The  importance  of  Rayleigh  scattering  and  the  surface  albedo  on 
the  solar  energy  absorbed  by  ozone  has  been  shown  by  Dave  and  Furukawa 
(1967).  Since  the  ozone  heating  rate  virtually  approaches  zero  at  about 
10  km,  ozone  was  modeled  as  an  absorbing  layer  overlaying  a  scattering 
layer.  In  a  clear  atmosphere,  Lacis  and  Hansen  (1974)  have  derived  a 
parameterized  equation  for  Rayleigh  scattering  albedo  at  the  top  of 
this  scattering  layer,  given  by 

=  0.28  /  (1  +  6.43  Pp)  .  (62) 


Using  this  expression,  the  total  reflection  at  the  top  of  this 
scattering  layer,  ^(Uq)*  including  the  surface  albedo,  may  be  deter¬ 
mined  by  a  ray  tracing  geometric  method.  As  shown  in  Fig.  8,  the  total 
reflection  at  the  top  of  the  Rayleigh  scattering  layer  is  given  by  an 
infinite  series  as  follows: 


-■(ko)  =  '"a  *  '•g  ‘a  "  ^<“0*  '"g  "'a  ‘a  "  ‘a<"o>  '"g  ^  ‘a 


+  t.(Wi)  '•f  • 


-  --(.(Mo)  +  [1  -  Tg  (1  -  r^)  [1  -  i-g  rf  +  .  .  .] 

=  ^(t'o'  *  tl  -  fjCuo)  ]  (1  -  f^)  /  (1  -  Tg  f^)  ,  (63) 

where  is  the  reflectivity  of  a  Lambert  reflecting  surface  and 


The  spherical  albedo  of  a  Rayleigh  scattering  atmosphere,  r  ,  is  defined 

G 

by 


|.i 

■0 


(65) 


The  numerical  value  of  r^  was  determined  to  be  0.0685. 

a 

Now,  the  rate  of  change  of  diffuse  upward  flux  at  any  z  >  z^  may 

a 

be  readily  calculated  by  replacing  the  zenith  angle  dependent  argu¬ 
ment  of  the  broadband  absorptivity,  (U3/PQ),  with  the  diffuse  path 
length,  (1.66  u^),  as  described  by  Goody  (1964).  Thus 

18  _ 

=  -  1.66  F  (z,)  I  f.  k.  exp  (-1.66  k.  u.)  ,  (66) 


where  the  upward  flux  at  z  =  z,  is  given  by 

0 

F”''(Za)  =  Sq  Uq  [1  -  A(u/Uq)]  r(uQ)  ,  (67) 

and  A(u/Mq)  is  the  total  gaseous  broadband  absorptivity  of  water  vapor, 
carbon  dioxide  and  ozone  in  the  column  from  the  top  of  the  atmosphere 
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The  instantaneous  ozone  heating  rate  at  any  height  above  the  top 


of  the  ozone  scattering  layer  is  given  by 


P  a 


M'S. 


P  a 


18 

r 

1=1 


1.66  F'^fz  )  y  f.  k.  exp  (-1.66  k.  n  M 
a  1  1 

(68) 

Ozone  heating  rates  were  calculated  for  a  clear  midlatitude  winter 
atmospheric  profile  (McClatchy  et  al.,  1972)  and  solar  zenith  angles  of 
0.2,  0.5,  and  1.0.  Figure  9  shows  the  results  obtained  from  the  param¬ 
eterization  model  in  the  dashed  curve  and  those  obtained  from  the  DOM 
in  the  solid  line  for  a  surface  albedo  of  zero.  The  two  methods  pro¬ 
duce  results  which  are  almost  indistinguishable  except  at  about  40  km 
where  the  parameterization  was  about  1°  C/24  hour  day  larger  than  the 
more  exact  method.  This  figure  illustrates  an  interesting  aspect  of 
ozone  heating  rates.  Even  though  the  maximum  ozone  concentration  is 
located  at  about  20-25  km,  the  maximum  heating  rate  occurs  at  about 
45-55  km,  depending  on  the  solar  zenith  angle.  As  Lacis  and  Hansen 
(1974)  have  explained,  this  is  caused  by  the  strong  absorption  of 
energy  in  the  ultraviolet  region  of  the  solar  spectrum  by  the  Hartley- 
Muggins  band  complex  at  altitudes  where  the  air  density  is  very  small. 

Results  in  Fig.  10  were  obtained  under  similar  circumstances  to 
those  of  Fig.  9,  except  the  surface  albedo  was  increased  to  100!^.  The 
sets  of  curves  are  again  in  remarkably  good  agreement.  It  should  he 
noteri,  howfver,  that  the  increased  surface  albedo  had  its  nwximum 
effect  in  the  10-35  km  range  where  the  diffuse  upward  radiation  con¬ 
tacted  the  maximum  ozone  concentration  in  the  atmosphere.  The  area  of 


maximum  ozone  heating  increases  only  slightly  by  increasing  the  surface 
albedo  from  zero  to  100%  because  the  Hartley-Huggins  bands  were  nearly 
saturated  by  the  absorption  of  the  direct  solar  flux. 


Hi 
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4.2.2  Water  Vapor  Absorptivity  and  Heating  Rate  Profiles 

The  water  vapor  broadband  absorptivity  is  also  described  by  Eq.  (59). 
Table  7  lists  the  wavelength  interval  over  which  Howard  et  al .  (1956) 
originally  measured  the  total  absorptivity  of  each  water  vapor  band. 

Also  listed  are  the  solar  flux  in  each  interval  after  Thekaekara  (1974) 
and  the  computed  fractional  solar  flux. 

In  Fig.  11,  the  broadband  absorptivity  for  water  vapor  path  lengths 
_2 

of  0.01  to  10  g  cm  ,  computed  at  STP  by  Eq.  (59),  are  compared  with 
other  published  data.  Two  important  points  are  present  in  this  figure. 
First,  the  broadband  absorptivity  curve  is  closest  to  the  line-by-line 
calculations  taken  from  Fouquart  and  Bonnel  (1980),  in  both  the  small 
and  large  water  vapor  path  lengths.  Secondly,  the  slopes  of  the  broad¬ 
band  and  line-by-line  absorptivity  curves  are  nearly  parallel.  Since 
the  heating  rate  is  directly  proportional  to  the  slope  of  the  absorp¬ 
tivity  curve,  the  heating  rate  curves  produced  by  the  broadband  water 
vapor  absorptivity  should  be  quite  accurate. 

The  rate  of  decrease  in  downward  flux  due  to  water  vapor  absorp¬ 
tivity  may  be  obtained  by  substituting  Eq.  (5)  into  Eq.  (55)  and  after 
differentiating  with  respect  to  water  vapor  path  length  u^,  we  obtain 


IL 

3Ui 


=  -^0  ''0  '“9io 


6  D.f . 
e  I 
i=l 


— 

+  1 

"l  -^i' 
.  p  ^ 
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(69) 


Table  7  Fractional  Solar  Flux 


1^0  band 

Interval  (iim) 

0.94 

0.8700-0.9900 

1.1 

1.0753-1.2048 

1.38 

1.2500-1.5385 

1.87 

1 .6949-2.0833 

2.70 

2.3041-2.9940 

3.20 

2.9940-3.5714 

in  H^O  Bands 

Solar  Flux  Fractional  Flux 
( W  m  ‘  ) 

102.810  0.0760 

71.488  0.0528 

99.034  0.0732 

52.256  0.0388 

32.733  0.0242 

11.921  0.0088 


wtiere  i  .  =  K./D^. .  Moreover,  from  Eq.  (8)  we  find 
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Broadband  water  vapor  absorptivity  compared 


exception  that  F  (z^)  is  replaced  by  F  (zq)  with  Zq  =  0.  It  must  also 
be  recognized  that  the  limits  of  integration  have  been  reversed  in  the 
reduced  pressure  expression.  The  total  heating  rate  at  level  z(Uj),  due 
to  water  vapor  absorption  of  direct  and  diffuse  solar  radiation  is  then 
given  by 


9t 


1 


[S  I  Y  (u./p  i  +  1.66  F^zJ  I  Y.  (1.66  u,)] 
^pf’a  ^  i=l  1  ^  ^  i=l  ^  ^ 


(73) 


Figure  12  depicts  the  exact  and  parameterized  water  vapor  heating 
rate  using  a  midlatitude  winter  profile  (McClatchey  et  al . ,  1972)  for  a 
clear  sky,  a  surface  albedo  of  zero,  and  solar  zenith  angles  of  0.2, 

0.5,  and  1.0  for  the  incident  solar  beam.  Setting  the  surface  albedo 
to  zero  produced  heating  rates  which  were  due  only  to  the  absorption 
of  direct  solar  flux.  In  general,  the  parameterization  and  the  exact 
curves  agree  within  about  0.1°  C  in  the  lower  troposphere  per  24  hrs 
of  heating.  Above  about  9  km,  the  disagreement  increases  to  a  maximum 
of  0.25°  C  per  24  hrs  with  the  sun  directly  overhead  for  the  entire 
time  period. 

In  Fig.  13,  the  surface  albedo  was  set  to  100%  to  maximize  the 
effect  of  absorption  of  the  diffuse  radiation  component  in  the  same 
atmospheric  profile.  A  comparison  of  this  curve  with  the  preceding 
one  reveals  the  effect  of  the  diffuse  component  was  confined  to  the 
lowest  6  km  of  the  troposphere  and  increased  toward  the  ground  where 
the  water  vapor  concentration  was  a  maximum.  Also,  the  magnitude  of 
the  heating  rate  caused  by  absorption  of  diffuse  radiation  was  consider¬ 
able,  approximately  tripling  at  the  surface  but  becoming  less  than  0.1°  C 


r  vapor  heating  rates  for  a  clear  midlatitude  winter  profile 
surface  albedo  r„  is  set  to  be  zero. 


at  6  km.  The  error  introduced  by  using  the  diffusivity  factor  to  repre¬ 
sent  the  scattering  over  all  zenith  angles  is  highlighted  by  the  param¬ 
eterized  heating  rates  becoming  larger  than  those  produced  by  the  exact 
method  in  the  lower  troposphere  when  the  surface  albedo  was  increased 
to  100%.  In  Fig.  12,  the  parameterized  method  produced  smaller  heating 
rates  throughout  the  same  region. 

Figure  14  compares  the  same  methods  for  a  clear  atmosphere  using  the 
midlatitude  summer  profile  compiled  by  McClatchey  et  al.  (1972)  with  a 
surface  albedo  of  15%.  This  albedo  is  an  average  albedo  for  the  earth, 
consequently  this  figure  shows  the  magnitude  of  heating  rate  errors 
which  will  typically  occur  using  the  parameterization  method. 


4.2.3  Carbon  Dioxide  Absorption  and  Heating  Rate  Profiles 

Since  the  spectral  interval  in  which  CO^  absorbs  solar  radiation 
in  the  2.7  pm  band  is  completely  contained  in  the  water  vapor  band,  the 
broadband  water  vapor  absorptivity  can  easily  be  modified  to  include 
its  effects.  As  previously  pointed  out,  the  contribution  of  carbon 
dioxide  to  the  total  atmospheric  absorptivity  in  the  2.7  pm  area  of  the 
electromagnetic  spectrum  increases  as  the  total  water  vapor  absorptivity 
decreases  in  this  area.  This  idea  is  represented  by 


A{H20  +  CO2)  = 


A^(H20)  dv  +  6 

J  V 


A^(C02)  dv  . 


V  =  2.7  pm 


(74) 


where  6  takes  a  value  between  0  and  1  as  depicted  in  Fig.  4. 

Using  the  concept  shown  in  Eq.  (59),  the  combined  water  vapor  and 
carbon  dioxide  broadband  absorptivity  is  then 


A(Ui2/Po)  ^  ^i  ^i^^^l^^O^  *  ^5  ^5 


(75) 
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midlatitude  summer  profile 


where  is  the  fractional  solar  flux  in  the  carbon  dioxide  3480-3800 
cm"^  frequency  interval  with  a  value  of  0.0073,  and  subscript  5  refers 
to  the  2.7  vim  band. 

Atmospheric  heating  caused  by  absorption  of  solar  flux  by  carbon 
dioxide  may  be  calculated  similarly  to  that  of  water  vapor,  using  Eqs. 
(57)  and  (58).  The  rate  of  change  of  downward  solar  flux  due  to  CO2 
absorption  is 


lEl  =  _i_ 

9U2  ~  3u2 


So  Mq  [1  -  6  fg  Ag  (u2/Pq)] 


■^0  ^^0 


6  Y(u2/Pq)  +  fg  Ag  (u2/^Jo^  9Uo 


(76) 


The  rate  of  change  of  weighting  factor  with  respect  to  carbon  dioxide 
path  length,  U2»  is  found  by  a  finite  difference  routine.  The  rate  of 
change  of  diffuse  upward  solar  flux  is  always  zero  since  the  total  water 
vapor  absorptivity  for  diffuse  radiation  in  the  2.7  urn  band  is  greater 
than  500  cm~^  so  6  is  always  zero.  This  occurs  because  the  water  vapor 
path  length  for  the  diffuse  beam  is  integrated  from  the  surface  to  the 
top  of  the  atmosphere.  Consequently,  the  CO2  heating  rate  equation  is 
given  by  substituting  Eq.  (76)  into  Eq.  (55): 


C  p  ^0  ^0  ^5  ^5 


(77) 


A  typical  heating  rate  curve  for  carbon  dioxide  is  shown  in  Fig. 

15  for  a  solar  zenith  angle  of  0.5.  Two  important  points  are  shown  in 
this  figure.  First,  there  is  no  heating  produced  below  about  12  km  where 
H2O  absorption  becomes  large  and  the  weight  factor  for  CO2  absorption 


approaches  zero.  Secondly,  the  maximum  carbon  dioxide  heating  occurs 
about  30  km  where  the  weight  factor,  6,  has  a  value  of  one  along  with 
small  air  densities. 

4.3  Total  Heating  Rate  and  Net  Flux  Profiles 

Climate  modeling  not  only  requires  accurate  heating  rate  computa¬ 
tion  but  also  accurate  net  fluxes,  especially  at  the  top  and  bottom  of 
the  atmospheric  model.  The  total  heating  rate  is  the  sum  of  the  indi¬ 
vidual  heating  rates  for  water  vapor,  ozone,  and  carbon  dioxide  at  each 
level.  However,  the  net  flux  is  complicated  by  Rayleigh  scattering. 
Molecular  scattering  tends  to  decrease  the  direct  solar  beam  and  in¬ 
creases  the  diffuse  upward  beam  as  can  be  seen  in  Fig.  16  which  depicts 
the  downward  and  upward  fluxes  for  both  methods  for  a  tropical  profile. 
The  zenith  angle  is  0.6  and  the  surface  albedo  is  0.9.  This  effect  is 
greatest  in  the  lower  troposphere  where  the  air  density  is  largest. 

It  is  apparent  that  the  net  fluxes,  defined  by  Eq.  (45),  for  these  two 
methods  disagree  by  about  50-75  W  m  throughout  the  profile,  even 
though  the  heating  profiles,  depicted  in  Fig.  17,  show  good  agreement. 
An  empirical  correction  must  be  made  to  the  net  flux  calculation  of  the 
parameterization  model  if  agreement  is  to  be  otained  with  a  more  exact 
model  which  includes  Rayleigh  scattering  effect. 

The  upward  flux  at  z,,  as  given  by  Eq.  (67),  is  denoted  here  as 

a 

F^(z-).  We  then  match  the  upward  flux  from  this  equation  with  the  up- 

a 

ward  flux  resulting  from  diffuse  absorption  by  water  vapor  in  the  area 
below  z..  The  latter  expression  is  given  by 

d 

=  F-‘(zJ  (1  -  -  Zj)  fg  (1  -  Jj)  .  (78) 
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The  empirical  correction  is  the  difference  between  these  two  fluxes, 
i  .e. , 

C  =  F^Zg)  -  F^z^)^  .  (79) 

The  upward  fluxes  at  all  heights  below  z^  are  now  adjusted  by  adding 

a 

this  correction  to  all  upward  fluxes  initially  calculated  using  Eq.  (78) 
The  result  of  computing  the  net  flux  profile  in  this  manner  is  demon¬ 
strated  in  Fig.  17  which  is  for  the  same  atmospheric  profile  used  in 
the  preceding  figure.  Good  agreement  is  now  observed  in  the  net  flux 
profile.  Extensive  experiments  utilizing  various  atmospheric  profiles, 
solar  zenith  angles  and  surface  albedos  also  show  good  agreement  (within 
about  5  W  m  )  for  net  flux  profiles  between  the  parameterization  and 


a  more  exact  method. 


Section  5 


THE  BROADBAND  APPROACH  IN  A 
CLOUDY  ATMOSPHERE 


5.1  General 

Most  parameterization  schemes  use  constants  or  some  form  of  the 
two-stream/Eddington  approximation  to  calculate  cloud  reflection  and 
transmission.  These  types  of  approximations  are  inherently  inaccurate 
but  their, use  is  justified  on  the  grounds  of  computational  simplicity. 
They  are  based  on  the  assumption  that  all  radiation  transmitted 
through  the  cloud  is  diffuse  and  fail  to  take  into  account  zenith 
angle  dependence  on  cloud  reflection  and  transmission. 

In  conjunction  with  gaseous  broadband  absorptivities  previously 
defined,  clouds  and  the  reflecting  surface  are  considered  as  a  system 
in  developing  parameterization  equations  for  fluxes  and  heating  rate 
calculations.  Direct  and  diffuse  broadband  solar  radiation  streams 
are  kept  separated  in  the  formulation  of  upward  and  downward  fluxes. 
Multiple  reflections  within  the  cloud-ground  system  are  handled  by  a 
modification  to  the  adding  method  using  diffuse  radiation  generation 
functions  which  are  defined  at  the  top  of  each  reflecting  layer  inside 
the  system.  The  heating  rate  within  a  cloud  layer  is  computed  using 
the  flux  density  divergence  for  the  layer  containing  the  cloud. 

5.2  Parameterization  of  Cloud  Optical  Properties 

Liou  and  Wittman  (1979)  have  shown  that  cloud  radiative  properties 
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(i.e.,  absorptivity,  transmissivity,  and  reflectivity)  are  adequately 
described  by  vertical  liquid  water  content  (W)  and  solar  zenith  angle 
(Mq)  once  the  cloud  microphysical  properties  have  been  established. 
However,  their  parameterization,  which  employed  16  coefficients  each 
for  cloud  reflection  and  transmission,  proved  to  be  cumbersome.  In 
addition,  their  definitions  were  not  in  complete  agreement  with  the 
broadband  gaseous  absorptivity  developed  in  this  study  as  pointed  out 
in  Section  3.  Thus,  a  new  cloud  parameterization  was  carried  out  for 
reflection  and  transmission.  After  trial  and  error  using  a  number  of 
known  mathematical  functions,  the  best  cloud  reflection  and  transmission 
parameterization  takes  the  form 


S(W,Mq)  = 


-1 

'o 


(SO) 


where 

4  „  1 

^(W)  =  I  >  n  =  1, 2,3,4  ,  (81) 

"  m=l 

and  S(W,^lQ)  refers  to  both  cloud  reflection  and  transmission.  In  these 
equations,  superscripts  n  and  m  on  )Jq  and  W  represent  the  index  power. 

This  parameterization  is  accomplished  in  the  following  manner. 
Cloud  optical  properties  are  obtained  for  five  cloud  thicknesses  and 
five  solar  zenith  angles  making  25  data  points  for  each  cloud  optical 
property.  The  number  of  data  points  needed  to  define  reflection  and 
transmission  is  reduced  significantly  from  the  90  used  by  Liou  and 
Wittman  (1979),  since  the  shape  of  the  curves  fitted  are  known  in 
advance.  The  cloud  radiative  property  is  plotted  against  the  solar 
zenith  angle  for  a  constant  vertical  liquid  water  content  as  depicted 
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in  Fig.  18  (the  stratus  cloud  transmission  is  employed  here  for  illus¬ 
tration).  This  produces  a  family  of  curves  which  is  then  fitted  by 
polynomials  described  in  Eq.  (80).  Next,  values  of  a^(W)  are  plotted 
against  W  as  shown  in  Fig.  19  (cumulus  reflection  is  used  for  illustra¬ 
tion)  and  this  family  of  curves  are  fit  by  a  polynomial  of  the  type 
given  in  Eq.  (81).  The  result  is  that  each  cloud  optical  property  is 
parameterized  using  six  to  nine  coefficients  except  for  cirrus  which 
require  16  coefficients.  The  accuracy  of  the  fit  for  each  cloud  radia¬ 
tive  property  is  shown  in  Table  8  along  with  the  cloud  thicknesses  and 
range  of  vertical  liquid  water  content  used  in  the  parameterization. 

This  parameterization  method  fits  the  cloud  reflection  generated  by 
the  discrete-ordinates  method  within  about  0.5%  on  the  average,  while 
the  cloud  transmission  is  fit  within  about  0.25%.  A  more  detailed 
comparison  of  the  cloud  radiative  properties  generated  by  the  discrete- 
ordinates  method  and  the  parameterization  are  given  in  Figs.  20  and  21. 
Figure  20  shows  the  comparison  between  the  altostratus  transmission  at 
various  W  and  while  Fig.  21  gives  the  comparison  between  the  cirrus 
reflection.  The  predictor  coefficients  given  by  Eq.  (81)  are  listed 
in  Tables  9-16  for  the  four  cloud  types. 

This  method  required  the  surface  albedo  to  be  set  to  zero  and 
assumed  that  Rayleigh  scattering  was  negligible,  which  it  is  for  all  clouds 
except  cirrus.  When  cirrus  is  the  only  cloud  present  in  the  atmos¬ 
phere  and  the  surface  albedo  is  zero,  the  effect  of  Rayleigh  scattering 
can  be  seen  by  summing  the  absorption,  transmission,  and  reflection 
for  each  case  as  shown  in  Table  17.  The  sum  should  be  one,  which  it 
is  for  all  other  cloud  types,  but  is  greater  than  one  in  all  cirrus 
cases.  This  occurs  because  in  clouds  with  high  transmission,  Rayleigh 
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Fig.  20  Comparison  of  the  exact  ( - )  and  parameterized  ( - )  methods 

for  the  altostratus  transmission  for  various  solar  zenith  angles 
and  liquid  water  contents. 


REFLECTION 


Table  9  Predictor  coefficients  for  the  stratus  reflection 


with  a  cloud  base  at  1  km. 


^n]m\^ 

1 

2 

3 

1.1770387 

-0.00117267 

-0.00155989 

^2,m 

-0.28671732 

1.3770553 

-0.08290901 

Table  10  Predictor  coefficients  for  the  stratus  transmission 
with  a  cloud  base  at  1  km. 


N.  m 

1 

2 

3 

•^l.m 

2.1467536 

3.0933553 

0.06150582 

^2,m 

-2.5417457 

-3.2657469 

-0.06695797 

^3,m 

1.2185517 

1.4555492 

0.03479150 

Table  11  Predictor  coefficients  for  the  c’jmulus  reflection 


with  a  cloud  base  at  2  km. 


Table  12  Predictor  coefficients  for  the  cumulus  trans¬ 
mission  with  a  cloud  base  at  2  km. 


b,  2.0945402  4.8603017  0.07852117 

1  ,m 

b„  „  -2.5646724  -5.332700  -0.11783940 

2,m 


b,  1.2151820  2.5195630  0.05836719 

3,m 


Table  13  Predictor  coefficients  for  the  altostratus 
reflection  with  a  cloud  base  at  4  kin. 


1.1818476  -0.069P8557  0.01?51b02 

0.74679249  4.7541147  -0.71275650 


Table  14  Predictor  coefficients  for  the  altostratus 
transmission  with  a  cloud  base  at  4  km. 


1 

2 

3 

^n^iri\^ 

^1  ,m 

2.1468193 

17.022049 

1.4969625 

*^2,m 

-2.8530089 

-18.961483 

-2.0434877 

1.4058295 

8.9499630 

1.0196345 

Table  15  Predictor  coefficients  for  the  cirrus  reflection  with  a 


cloud  base  at  8  km. 


m 

1 

2 

3 

4 

b 

n,m 

,m 

0.875416 

0.856228x10^ 

-0.361108x10^ 

0.393977x10^ 

^2,m 

0.185111x10^ 

-0.855647x10^ 

0.354341x10^ 

-0.386822x10^ 

^3,m 

0.667364x10^ 

0.117047x10^ 

-0.595519x10^ 

0.689059x10^ 

^4  ,m 

-0.622586x10^ 

-0.529730x10^ 

0.313787x10^ 

-0.376506x10^ 

Table  16  Predictor  coefficients  for  the  cirrus  transmission  with 
a  cloud  base  at  8  km. 
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scattering  will  reflect  flux  arriving  below  the  cloud  upward  again. 

In  the  discrete-ordinates  method,  this  additional  upward  flux  cannot  be  dis¬ 
tinguished  from  flux  originally  transmitting  through  the  cloud.  However, 
f'l.  Rayleigh  ‘^';atter;ng  on  cima-  c’oul':.  .  .  ;  i/ed 

in  the  same  way  that  cloud  reflection  and  transmission  are  parameterized 
for  all  cloud  types.  Figure  22  depicts  the  inverse  of  the  Rayleigh 
effect  on  cirrus  clouds  plotted  against  the  solar  zenith  angle  for  a  con¬ 
stant  liquid  water  content.  Coefficients  ap(W)  in  Eq.  (80)  are  again 
obtained  by  fitting  this  family  of  curves  by  a  third-order  polynomial 
in  the  zenith  angle.  Figure  23  shows  the  plot  of  these  coefficients 

versus  liquid  water  content  from  which  b  in  Eq.  (81)  are  obtained 

n  ,m 

after  fitting  this  family  of  curves  with  the  same  type  of  polynomials. 

The  accuracy  of  this  fit  on  the  average  is  within  0.026%  of  the  values 
obtained  by  the  discrete-ordinates  method.  The  coefficients,  b 

n,m 

are  listed  in  Table  18. 

Several  points  should  be  understood  when  applying  these  predic¬ 
tion  coefficients  listed  in  Tables  9  through  17.  First,  they  apply 
to  solar  zenith  angles  within  the  range  0.2  -  1.0.  Secondly,  for  the 
reflectic-  of  stratus,  cumulus,  and  altostratus  clouds,  they  predict 
b2^  not  b2.  It  was  necessary  to  predict  b^^  if  the  number  of  coeffi¬ 
cients  used  to  predict  reflection  is  to  be  kept  to  six. 

An  inflection  point  exists  at  Uq  =  .2  for  most  cloud  types;  thus, 
large  errors  could  result  if  these  coefficients  are  applied  with  solar 
zenith  angles  less  than  0.2.  Errors  are  also  probable  when  used  wHh 
cloud  thicknesses  smaller  than  those  shown  in  Table  17  but  only  small 
errors  would  result  by  using  thicker  clouds.  The  vertical  liquid 
water/ice  content  of  the  cloud  is  given  by 
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SOLAR  ZENITH  ANGLE 


Fig.  22  Inverse  of  the  Rayleigh  scattering  effect  on  cirrus  clouds 

versus  solar  zenith  angle  for  a  constant  liquid  water  content 


Table  18  Predictor  coefficients  due  to  the  Rayleigh  effect  with  a 


cloud  base  at  8  km. 


'  1 

N.  m 

1 

2 

3 

4 

0.543956x10^ 

0.330860x10^ 

-0.705047x10^ 

0.625258x10® 

^2,m 

0.500547x10^ 

. 

-0.126917x10^ 

0.305654x10^ 

-0.274697x10^ 

^3,m 

-0.492899x10^ 

0.138483x10^ 

-0.349088x10^ 

0.321333x10^ 

^4,m 

0.231012x10^ 

-0.422886x10^ 

0.104622x10® 

-0.100449x10^ 

W  =  w  Az  , 


(82) 


where  Az  is  the  cloud  thickness  in  meters  and  w  is  0.78,  0.33,  0.24, 

_3 

and  0.00336  g  m  for  stratus,  cumulus,  altostratus,  arc*  '■irrus, 

2  -2 

respectively.  The  units  of  M  are  in  10  g  m  . 

Additionally,  a  temperature  correction  to  the  cloud  reflection 
and  transmission  can  be  made  using  the  form, 

S(w,Pq,T)  =  S(w,Pq)  +  C(T)  ,  (83) 

where  S(w,Pq)  is  determined  from  Eq.  (80)  and  C(T)  the  temperature  cor¬ 
rection  depends  on  the  mean  temperature  of  the  cloud.  The  temperature 
correction  to  cloud  reflection  and  transmission  is  an  important  new 
concept  in  the  parameterization  of  cloud  radiative  properties  because 
it  allows  the  mean  temperature  and  consequently  the  base  of  the  cloud 
to  vary. 

The  correction,  C(T),  was  found  for  a  40  degree  temperature  range 
for  each  cloud  type  considered  in  the  research.  The  temperature  correc¬ 
tion  was  made  by  varying  the  temperature  of  the  cloud  in  the  more  exact 
method  and  plotting  the  resulting  radiative  property  against  the  mean 
cloud  temperature  for  a  given  solar  zenith  angle  as  shown  in  Fig.  24. 
Since  each  curve  closely  approximates  a  straight  line,  the  parameter¬ 
ization  could  easily  be  made  by  describing  the  slope  of  each  curve  as 
a  function  of  zenith  angle  for  a  given  vertical  liquid  water  content, 
as  shown  in  Fig.  25.  The  result  is  a  temperature  coefficient  of  the 
form 

3  _  , 

C(w,u^,T)  =  (T  -  T^g^)  I  c^(W)  ,  (84) 

n=l 
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where 


Cn(W) 


3 

I  d 

m=l 


,m-l 


n,m 


(85) 


T  is  the  mean  temperature  of  the  cloud,  and  is  determined  from 
the  reference  cloud  base  temperature  for  each  cloud  type.  These  refer¬ 
ence  temperatures  are  265.2,  255.7,  231.7,  and  268.7°  K  for  cumulus, 
altostratus,  cirrus,  and  stratus,  respectively.  These  predictor  co¬ 
efficients  and  applicable  temperature  ranges  are  listed  in  Tables  19  - 
26.  The  result  of  the  cloud  radiative  parameterization  is  that  once 
the  cloud  type,  thickness,  and  base  are  fixed,  the  reflection  and  trans¬ 
mission  are  defined  for  given  zenith  angles.  Deviations  of  the  param¬ 
eterized  values  from  those  of  the  more  exact  method  when  the  tempera¬ 
ture  correction  is  made  are  within  the  range  shown  in  Table  8. 


5.3  Heating  Rates  and  Net  Fluxes  in  Cloudy  Atmospheres 

Heating  rate  calculations  for  an  atmosphere  containing  a  single 
cloud  are  divided  into  three  general  areas:  the  area  above  the  cloud 
top,  the  area  below  the  cloud  base,  and  the  cloud  layer  as  shown  in 
Fig.  26.  In  the  formulation,  direct  and  diffuse  components  of  the 
solar  flux  contribute  to  the  net  flux  and  heating  rate  at  each  level  in 
the  cloud  layer. 

The  cloud  has  a  transmission,  tj(pQ)  and  reflection,  which 

are  functions  of  the  vertical  liquid  water  content  and  solar  zenith 
angle  obtained  for  an  identical  cloud  when  the  surface  albedo  is  set 
to  zero.  Additionally,  the  cloud  has  a  global  (or  average)  reflection 
r^,  and  global  transmission,  t^,  with  respect  to  diffuse  radiation. 
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Table  19  Predictor  coefficients  for  the  temperature  correction  to 

the  cumulus  reflection.  The  temperature  range  used  is  from 
250  to  290“  K. 


*n,nNv 

1 

2 

3 

r 

\ 

J, 

1  ,m 

1.66157x.O"^ 

-2.87867x10"^ 

1.55598x10"^ 

-4.25249x10"^ 

-3.18357x10"^ 

2.83799x10"^ 

3.97082x10"^ 

6.38750x10"^ 

7.32911x10"^ 

-  ’ 

Table  20  Predictor  coefficients  for  the  temperature  correction  to 

the  cumulus  transmission.  The  temperature  range  used  is  from 
250  to  290“  K. 


Table  21  Predictor  coefficients  for  the  temperature  correction  to  the 
altostratus  reflection.  The  temperature  range  used  is  from 
245  to  285“  K. 


N.  m 

n,m\. 

1 

2 

3 

5.80363x10'^ 

2.00546x10'^ 

-5.73065x10'^ 

‘^2,m 

-6.03168x10'^ 

-2.35728x10'^ 

4.12499x10'^ 

^3,m 

3.57757x10'^ 

9.05252x10'^ 

-1.28976x10'^ 

Table  22  Predictor  coefficients  for  the  temperature  correction  to  the 
altostratus  transmission.  The  temperature  range  used  is  from 
245  to  285“  K. 


m 

1 

2 

3 

^l,m 

-6.24068x10'® 

-1.93991x10'® 

5.20037x10'® 

‘^2,m 

-3.48751x10'^ 

2.56122x10'^ 

-4.61342x10'® 

‘*3,m 

1.02351x10'^ 

-7.42415x10'® 

1.32371x10'® 

Table  23  Predictor  coefficients  for  the  temperature  correction  to  the 
cirrus  reflection.  The  temperature  range  used  is  from  205 
to  245®  K. 


-4.80458x10' 

2.39591x10' 


1.84610x10 


-5.23109x10 

7.74293x10' 


2.96808x10 


4.93013x10 


-4.28070x10 


2.89473x10 


Table  24  Predictor  coefficients  for  the  temperature  correction  to  the 
cirrus  transmission.  The  tempe»‘ature  range  used  is  from  205 
to  245®  K. 


N.  m 

1 

2 

3 

-1.69587x10'^ 

4.45949x10'^ 

-4.90636x10'^ 

^2,m 

-4.45260x10'^ 

-3.44686x10'^ 

4.84774x10'^ 

*’3,m 

1.71100x10'^ 

1.75961x10'^ 

-3.05332x10'^ 

Table  25  Predictor  coefficients  for  the  temperature  correction  to  the 
stratus  reflection.  The  temperature  range  used  is  from  260 
to  300*  K. 


N.  m 

1 

2 

3 

-5.50859x10"^ 

^l.m 

1.97809x10"^ 

-4.27582x10"^ 

^2.m 

-3.35439x10"^ 

-1.52791x10"^ 

1.34508x10"^ 

^3,m 

2.44638x10"'^ 

4.42694x10"^ 

-3.42760x10"^ 

Table  26  Predictor  coefficients  for  the  temperature  correction  to  the 
stratus  transmission.  The  temperature  range  used  is  from  260 
to  300°  K. 


-1.11568x10 


The  Lambert  reflecting  surface  at  z  =  0,  has  a  value  of  r^  and  subscripts 
t  and  b  represent  the  cloud  top  and  base,  respectively. 

Let  aj^  be  the  total  gaseous  absorptivity  between  the  reflecting 
surface  and  cloud  base,  and  a^  the  diffuse  absorptivity  in  the  same 
layer.  Define  62  as  the  "generation  function"  which  generates  the  total 
upward  diffuse  flux  at  the  surface  and  assign  the  downward  flux  at  the 
cloud  top  a  value  of  unity  for  simplicity,  i.e.,  F^(z^)  =  1.  Let  the 
upward  flux  at  z  =  0  after  the  first  reflection  be  D.  Then,  the  total 
upward  flux  at  z  =  0,  defined  by  G2,  is  given  by  the  infinite  series, 

G2  =  D  +  D  Tg  7*  (1  -  1^)2  +  D  r^  7*^  (1  -  +  .  .  .  (86) 

where  D,  as  shown  in  Fig.  27  is  given  by 

D  =  tj(yQ)  (1  -  aj)  rg  .  (87) 

The  star  in  Eq.  (86)  denotes  reflected  radiation  originating  below  the 
cloud.  This  infinite  series  converges  to 

62  =  D  r^  7*  (1  -  aj)^  [1  -  (1  -  a^)^  Cg  r*]'^  .  (88) 

If  $2  is  defined  as 

$2  =  (1  -  a^)^  rg  7*  [1  -  (1  -  a^)^  Cg  r*]'^  ,  (89) 

then  G2  can  be  expressed  by 

62  =  tj(yQ)  (1  -  a^)  r  (1  +  $2) 

^  »  (90) 

A  B 

where  term  A  represents  the  direct  solar  flux  reaching  the  reflecting 


surface  and  term  B  contains  the  multiple  reflection  effects  which  occur 
between  the  cloud  and  the  reflecting  surface.  Since  the  downward  flux 
at  the  cloud  top  is  assumed  to  be  unity,  the  upward  flux  at  the  surface 

F^Zq)  =  62- 

The  upward  flux  at  the  cloud  base  is  the  total  upward  flux  at  the 
earth's  surface  attenuated  by  the  diffuse  absorptivity  of  the  atmospheric 
layer  between  Zq  and  the  cloud  base.  Consequently, 

F^z^)  =  G2  (1  -  ij)  .  (91) 

The  downward  flux  arriving  at  the  cloud  base  has  two  components.  The 
first  is  direct  solar  flux  transmitted  through  the  cloud  and  the  second 
represents  diffuse  solar  flux  originating  at  Zq  and  reflected  by  the 
cloud.  This  solar  flux  is  represented  by 

F^Zjj)  =  tj(yg)  +  G2  (1  -  aj)  ?*  .  (92) 

The  total  downward  iux  at  Zq  is  Eq.  (92)  with  each  component  atten¬ 
uated  by  absorption  in  the  layer  between  Zq  and  the  cloud  base.  Thus, 

F'^(Zq)  =  tj(pQ)  (1  -  a^)  +  G2  (1  -  a^)^  .  (93) 

The  total  upward  flux  at  the  cloud  top  is  the  sum  of  the  direct  down¬ 
ward  flux  reflected  by  the  cloud  and  additional  diffuse  solar  flux 
which  originates  at  Zq  and  is  subsequently  transmitted  through  the 
cloud  during  the  multiple  reflection  process.  This  final  expression 
is  given  by 

F^z^)  =  r*(MQ)  =  +  G2  (1  -  a^)  t*  .  (94) 
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The  heating  rate  at  the  middle  of  the  cloud  is 


£[  -  _1_  AE 
at  ~  ^  -  Az  ’ 


where  the  radiative  flux  density  divergence  is  given  by 


AF  -  [F^z^)  -  F^z^]  -  [F^zjj)  -  F^Zj^)]  , 


(95) 


(96) 


and  p,  is  the  mean  air  density  of  the  cloud  layer  and  Az  is  the  cloud 

a 

thickness. 

Since  the  upward  and  downward  fluxes  at  the  cloud  base  contain 
direct  and  diffuse  components,  heating  rates  in  the  clear  atmosphere 
are  determined  as  previously  described  for  a  cloud  free  atmosphere 
with  the  exception  that  the  upward  flux  at  level  z  above  the  cloud 
top  is 

F^(z)  =  Sq  Pq  [1  -  A(z^)]  r*(pQ)  [1  -  A(z)] 

z  >  z^  (97) 

where  A(z)  =  A  (1.66  z).  On  the  other  hand,  the  downward  and  upward 
fluxes  at  level  z  below  the  cloud  base  takes  the  form 

F’‘(z)  =  Sq  Pq  [1  -  A(z^)]  {  tj(PQ)  [1  -  A{z)] 

^  (1  -  a^)  r*  [1  -  J{z^  -  z)] 

F^z)  =  Sq  Pq  [1  -  A(z^)]  (1  -  a^)  [1  -  A(z)]  . 

z  Wb 


Figures  28,  29  and  30  show  heating  rates  for  (.unuK, 
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and  cirrus  clouds,  respectively  In  a  mldlatltude  winter  atmosphere. 

Each  cloud  Is  one  km  thick  and  the  sun  directly  overhead.  The  surface 
albedo  was  set  to  zero  In  all  cases  to  maximize  the  effect  of  the 
absorption  of  the  direct  solar  flux  on  heating  rates  In  the  three  areas 
of  the  model.  Important  cloud  properties  are  listed  In  Table  27. 

Heating  rates  above  and  below  the  cloud  are  generally  quite  good 
for  all  three  cloud  cases.  Poorer  agreement  Is  found  Immediately  above 
the  altostratus  and  cumulus  clouds  where  the  diffuse  component  dominates 
the  heating.  The  reflection  values  of  these  clouds  when  Mq  *  1 
0.823  and  0.763,  respectively.  Below  the  cirrus  cloud,  the  parameter¬ 
ized  model  produces  heating  rates  which  are  generally  larger  than  the 
discrete-ordinates  model.  This  results  from  a  more  accurate  treatment 
of  Rayleigh  scatter  In  the  exact  method.  Rayleigh  scattering  reduces 
the  downward  flux  at  each  level  below  the  cloud  as  seen  In  Fig.  16, 
thus  reducing  the  heating  rate  In  each  layer. 

The  effect  of  Rayleigh  scatter  was  also  analyzed  In  each  cloud 
layer  and  found  to  be  most  Important  for  clouds  which  are  semi-transparent 
such  as  cirrus.  This  Information  Is  also  shown  In  Table  27.  The  up¬ 
ward  flux  densities  at  the  base  of  the  cumulus,  altostratus,  and  cirrus 
clouds  were  5.3,  3.0  and  42.5  W  m  ,  respectively.  In  each  of  the  cases. 
The  comparable  component  In  the  parameterization  model  was  zero  In  each 
case.  Since,  In  the  more  exact  model,  energy  absorbed  by  the  cumulus 
cloud  was  136  W  m”^,  164  W  m~^  for  altostratus,  and  83  W  m”^  for  the 
cirrus  cloud.  It  Is  obvious  that  Rayleigh  scatter  Is  very  Important  in 
cirrus  clouds  while  Insignificant  In  the  cumulus  and  altostratus  clouds. 

If  the  Rayleigh  component  Is  Ignored,  the  cloud  heating  rates 
agree  within  0.36“  C/24  hr  for  cumulus  and  0.20“  C/24  hr  for  altostratus. 
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However,  the  cirrus  heating  rate  is  50%  lower  in  the  parameterization 
model.  Adding  the  Rayleigh  component  to  the  parameterized  model  results 
in  exact  agreement  in  the  first  two  cloud  types  and  only  a  0.2**  C/24 
hr  error  in  the  cirrus  cloud  case.  Consequently,  this  Rayleigh  scat¬ 
tering  correction  was  parameterized  using  Eq.  (80)  for  cirrus  clouds 
and  the  effect  is  added  to  the  absorptivity  of  the  cirrus  cloud. 

Increasing  the  surface  albedo  allows  the  multiple  reflections  and 
diffuse  radiation  effects  on  cloud  heating  rates  and  net  flux  profiles 
to  be  studied.  In  Figs.  31  through  33,  heating  rate  and  net  flux 
profiles  are  shown  for  a  constant  surface  albedo  of  15%  (global  albedo). 
Table  28  lists  heating  rates  for  the  exact  and  parameterized  models  as 
a  function  of  surface  albedo.  In  all  cases,  the  solar  zenith  angle  is 
one  and  each  cloud  is  1  km  thick.  Agreements  of  heating  rate  and  net 
flux  profiles  between  the  parameterization  method  and  a  more  exact 
method  are  generally  excellent.  There  is  an  area  of  disagreement 
between  the  two  models  in  the  8-12  km  region  of  the  heating  rate  curves 
where  the  more  exact  model  produces  slightly  larger  results.  Cumulus 
and  altostratus  net  flux  profiles  contain  no  flux  correction;  but, 
as  in  the  clear  atmosphere,  Rayleigh  scattering  in  a  cirrus  cloudy 
atmosphere  requires  a  correction  to  generate  accurate  net  flux  profiles. 
In  the  cirrus  case,  the  upward  Rayleigh  component  is  calculated  using 
Eq.  (80)  and  the  upward  flux  at  each  level  between  the  cloud  base  and 
the  surface  is  inferred  as  is  done  in  the  clear  sky  case. 

Analysis  of  the  data  presented  in  Table  28  reveals  that  this 
method  of  accounting  for  multiple  reflections  between  the  cloud  and 
ground  actually  overestimates  the  effects  which  occur  in  the  discrete- 
ordinates  method.  The  difference  between  the  heating  rates  produced 


32  Heating  rate  and  net  flux  profiles  in  a  midlatitude  winter  atmosphere  contain 
ing  a  single  altostratus  (4-5  km)  cloud. 


and  net  flux  profiles  in  a  midlatitude  winter  atmosphere  cont 
cirrus  (8-9  km)  cloud. 


by  the  exact  and  parameterization  models  for  cumulus  and  altostratus 
clouds  when  the  surface  albedo  was  zero  has  been  reduced  by  increasing 
the  ground  albedo  to  15%.  If  the  albedo  is  Increased  to  80%  in  the 
altostratus  case,  the  parameterization  model  produces  a  heating  rate 
about  one  degree  larger  than  the  exact  model.  The  error  in  the  multiple 
reflection  effect  is  more  clearly  seen  by  observing  this  difference  in 
the  cirrus  cloud  case.  The  difference  grows  from  -0.2,  to  -2.85°  C/24 
hr  as  the  ground  albedo  is  increased  from  zero  to  100%.  The  difference 
in  the  multiple  reflection  effect  of  the  parameterization  model  probab' 
results  from  a  complicated  relationship  between  ground  albedo,  atmos¬ 
pheric  absorptivity  of  diffuse  radiation,  and  Rayleigh  scattering. 


Section  6 


SUMMARY  AND  CONCLUSIONS 

It  was  previously  stated  that  the  parameterization  approach  to 
solving  the  solar  radiative  transfer  equation  is  necessary  because  the 
more  exact  approximations,  such  as  the  discrete-ordinates  method,  are 
far  too  time  consuming  and  expensive  to  use  in  global  circulation 
models  of  the  earth-atmosphere  system.  However,  many  of  the  earlier 
attempts  at  parameterizing  the  solar  radiative  transfer  equations  were 
too  simple  and  known  relationships,  such  as  solar  zenith  angle  and 
liquid  water  content  dependence  on  the  reflection  and  transmission 
of  clouds,  were  not  accurately  considered.  The  primary  objective  of 
this  reasearch  was  to  develop  a  parameterization  model  whit.i  would 
generate  reliable  and  economical  solar  heating  rate  and  net  flux 
profiles  for  inclusion  into  a  global  circulation  model.  The  param¬ 
eterization  model  developed  in  this  report  is  based  on  physical  and 
mathematical  principles  which  govern  radiative  transfer  in  planetary 
atmospheres  under  clear  and  cloudy  sky  conditions. 

First,  consider  the  theoretical  principles  upon  which  this  param¬ 
eterization  model  is  based.  Under  clear  sky  conditions,  water  vapor 
heating  rates  are  assumed  to  result  from  gaseous  absorption  only. 

This  assumption  is  consistent  with  other  parameterization  models. 
Heating  rate  comparisons  between  the  parameterization  and  discrete- 
ordinates  methods  generally  agreed  within  0.2®  C/24  hr  for  all  solar 
zenith  angles  and  any  surface  albedo.  Ozone  heating  rates  are  known 


to  be  dependent  upon  Rayleigh  scatter  and  the  surface  albedo.  Conse¬ 
quently,  ozone  was  modeled  as  an  absorbing  layer  overlaying  a  reflect¬ 
ing  layer  as  described  by  Dave  and  Furukawa  (1967),  Lacis  and  Hansen 
(1974),  and  others.  Ozone  heating  rate  comparisons  between  the  more 
exact  and  parameterization  models  were  very  good  with  about  1°  C/24  hr 
being  the  largest  disagreement  in  a  profile  which  had  a  maximum  heating 
of  35“  C/24  hrs.  It  was  found  that  a  correction  to  the  net  fluxes  in 
a  clear  atmosphere  must  be  made  if  profiles  generated  by  the  two  models 
are  to  agree  well.  The  correction  is  necessary  because  Rayleigh  scat¬ 
ter  depletes  the  downward  direct  solar  beam  and  increases  the  upward 
diffuse  solar  beam  in  the  troposphere  where  molecular  number  densities 
are  large. 

In  cloudy  atmospheres,  this  parameterization  model  allows  much 
more  flexibility  to  designers  of  general  circulation  models  while 
maintaining  its  strong  theoretical  basis.  Most  parameterization  models 
assign  constants  to  cloud  reflections  and  transmissions  or  use  a  form 
of  the  two-stream  approximation  to  calculate  these  properties  from  the 
optical  depth  of  the  cloud.  Although  the  latter  approximation  is  much 
more  advanced  than  assigning  climatological  values  to  clouds  based  upon 
their  type,  it  does  not  allow  the  solar  zenith  angle  dependence  upon 
cloud  reflection  and  transmission,  which  can  be  rather  large,  to  be 
taken  into  account.  This  parameterization  method  does  consider  this 
dependence.  Additionally,  cloud  radiation  properties  are  parameterized 
as  a  function  of  liquid  water/ice  content.  Since  global  circulation 
models  now  forecast  relative  humidity  and  temperature  fields,  it  is 
highly  desirable  that  the  parameterization  of  cloud  radiative  properties 
also  possess  a  temperature  dependence.  This  parameterization  is  the 
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only  model  known  of  at  this  time  which  accounts  for  the  temperature 
dependence  on  the  cloud  radiative  properties. 

The  modified  adding  method  used  in  this  parameterization  model 
treats  direct  and  diffuse  radiative  streams  separately  which  allows 
highly  accurate  heating  rates  to  be  computea  outside  the  clouds. 

Multiple  reflection  inside  the  cloud-reflection  surface  system  have 
been  considered  and  internal  net  fluxes  are  calculated  at  each  level 
within  the  system.  It  was  noted  that  the  adding  method,  in  its  present 
form,  tends  to  overestimate  the  multiple  reflections  effect  on  the 
cloud  heating  rate  in  an  atmosphere  containing  a  single  cloud.  However, 
the  error  is  not  significant  when  dealing  with  normal  surface  albedos 
in  the  0.07-0.5  range. 

Next,  consider  the  computational  efficiency  of  the  parameteriza¬ 
tion  model  against  that  of  the  discrete-ordinates  method.  Table  29 
depicts  computational  statistics  for  the  discrete-ordinates  method 
and  the  parameterization  method  developed  in  this  research  project. 

The  exact  method  used  42  spectral  band  intervals  to  describe  the  solar 
spectrum.  It  was  run  for  one  solar  zenith  angle  and  for  various  numbers 
of  layers  in  the  atmospheric  profile.  One  cloud  was  present  in  all 
the  models  shown  in  the  table.  The  results  readily  demonstrate  the 
need  for  parameterizing  the  radiative  processes  in  GCM  where  solar 
heating  rate  calculations  are  needed  at  each  grid  point  and  for  each 
time  step.  As  the  number  of  layers  in  the  atmospheric  profile  in¬ 
creases  the  core  size  and  CPU  time  required  to  compute  the  solar 
heating  increases  dramatically  as  in  a  non-linear  fashion.  The  exact 
model  will  compute  the  solar  heating  rate  in  10  layers  for  a  single 
grid  point  and  a  single  zenith  angle  in  approximately  five  minutes 
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using  nearly  60  K  of  storage  in  the  UNIVAC  1110  here  at  the  University 
of  Utah.  Increasing  to  15  layers  almost  triples  the  CPU  time  required 
for  the  computations.  This  occurs  because  the  unknowns  in  the  system 
of  equations,  described  in  Section  three,  increase  from  40  to  60.  The 
core  size  increases  significantly  since  numerous  arrays  are  also  increased 
from  40  to  60.  For  exactly  the  same  reasons,  increasing  the  number  of 
layers  in  the  atmosphere  to  25  layers  requires  nearly  an  hour  of  CPU 
time  at  well  over  100  K  to  solve  the  system  of  equations  for  only  one 
solar  zenith  angle.  Because  a  minimum  of  20-25  layers  are  needed  to 
adequately  define  the  solar  heating  rate  up  to  60  km  in  a  cloudy  atmos¬ 
phere,  the  exact  method  obviously  is  of  little  value  in  a  global  model. 

Since  the  parameterization  model  uses  analytical  functions  rather 
than  the  computationally  slower  and  more  inaccurate  finite  difference 
schemes,  it  can  compute  solar  heating  rates  in  a  60  layered  atmospheric 
model  in  only  about  nine  CPU  seconds  at  less  than  40  K  on  the  UNIVAC 
1110.  With  a  nearly  linear  relationship  between  the  number  of  calcu¬ 
lations  and  the  number  of  layers  in  the  parameterization  model,  it  is 
estimated  that  the  parameterization  model  would  execute  in  about  3  CPU 
seconds  if  the  number  of  layers  were  reduced  to  20.  Similarly,  the 
CPU  should  drop  to  about  13  K  storage.  Clearly,  the  parameterization 
model  computes  solar  heating  rates  several  orders  of  magnitude  faster 
than  the  exact  method. 

The  parameterization  model  is  soundly  based  upon  radiative  prin¬ 
ciples  which  govern  planetary  atmospheres  and  is  considerably  more 
physically  reliable  than  parameter! zations  currently  being  used  in 
global  circulation  models  especially  when  clouds  are  involved.  The 
basic  contribution  of  this  research  work  is  that  maximum  flexibility 
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1s  built  into  the  cloud  radiative  properties  since  both  cloud  transmis 
Sion  and  reflection  remain  functions  of  three  independent  cloud 
parameters:  vertical  liquid  water/ice  content,  solar  zenith  angle, 
and  mean  cloud  temperature.  Heating  rates  are  computed  by  analytical 
expressions  at  precise  points  rather  than  by  finite  difference  approxi¬ 
mations.  Finally,  separating  the  direct  and  diffuse  solar  flux  allows 
accurate  heating  rates  to  be  computed  above,  below,  and  between  clouds 
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